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Abstract 

Purpose - This paper aims at assessing the effect of (1) the statical 
admissibility of the recovered solution; (2) the ability of the recovered so- 
lution to represent the singular solution; on the accuracy, local and global 
cffcctivity of recovery-based error estimators for enriched finite element 
methods (e.g. the extended finite element method, XFEM). 
Design/methodology/approach - We study the performance of two 
recovery techniques. The first is a recently developed superconvergent 
patch recovery procedure with equilibration and enrichment (SPR-CX). 
The second is known as the extended moving least squares recovery (XMLS), 
which enriches the recovered solutions but does not enforce equilibrium 
constraints. Both are extended recovery techniques as the polynomial ba^ 
sis used in the recovery process is enriched with singular terms for a better 
description of the singular nature of the solution. 

Findings - Numerical results comparing the convergence and the effectiv- 
ity index of both techniques with those obtained without the enrichment 
enhancement clearly show the need for the use of extended recovery tech- 
niques in Zienkiewicz-Zhu type error estimators for this class of problems. 
The results also reveal significant improvements in the effectivities yielded 
by statically admissible recovered solutions. 

Originality /value - This work shows that both extended recovery pro- 
cedures and statical admissibility are key to an accurate assessment of the 
quality of enriched finite element approximations. 

Keywords extended finite element method; error estimation; linear elastic frac- 
ture mechanics; statical admissibility; extended recovery 
Paper Type Research Paper 

1 Introduction 

Engineering structures, in particular in aerospace engineering, are intended to 
operate with flawless components, especially for safety critical parts. However, 
there is always a possibility that cracking will occur during operation, risking 
catastrophic failure and associated casualties. 

The mission of Damage Tolerance Assessment (DTA) is to assess the influ- 
ence of these defects, cracks and damage on the ability of a structure to perform 
safely and reliably during its service life. An important goal of DTA is to esti- 
mate the fatigue life of a structure, i.e. the time during which it remains safe 
given a pre-existing flaw. 



Damage tolerance assessment relies on the ability to accurately predict crack 
paths and growth rates in complex structures. Since the simulation of three- 
dimensional crack growth is either not supported by commercial software, or 
requires significant effort and time for the analysts and is generally not cou- 
pled to robust error indicators, reliable, quality-controlled, industrial damage 
tolerance assessment is still a major challenge in engineering practice. 



The extended finite element method (XFEM) (Moes et al. 1999 1 is now one 



of many successful numerical methods to solve fracture mechanics problems. 
The particular advantage of the XFEM, which relies on the partition of unity 
(PU) property (Melenk and Babuska 19961 of finite element shape functions, 
is its ability to model cracks without the mesh conforming to their geometry. 
This allows the crack to split the background mesh arbitrarily, thereby leading 
to significantly increased freedom in simulating crack growth. 

This feat is achieved by adding new degrees of freedom to (i) describe the 
discontinuity of the displacement field across the crack faces, within a given 
element, and (ii) reproduce the asymptotic fields around the crack tip. Thanks 
to the advances made in the XFEM during the last years, the method is now 
considered to be a robust and accurate means of analysing fracture problems. 



has been implemented in commercial codes ( ABAQUS 2009 CENAERO 2011 1 



and is industrially in use for damage tolerance assessment of complex three 



dimensional structures ( Bordas and Moran 20061 Wyart et al. 2007 Dufiot 



2007). 



While XFEM allows to model cracks without (re)meshing the crack faces 
as they evolve and yields exceptionally accurate stress intensity factors (SIF) 
for 2D problems, [Bordas and Moran] ( |2006[ ), [Wyart et alT] ( |2007[ ), [Dufiot [ ( [2007[ ) 
show that for realistic 3D structures, a "very fine" mesh is required to accurately 
capture the complex three-dimensional stress field and obtain satisfactory stress 
intensity factors (SIFs), the main drivers of linear elastic crack propagation. 

Consequently, based on the mesh used for the stress analysis, a new mesh 
offering sufficient refinement throughout the whole potential path of the crack 
must be constructed a priori, i.e. before the crack path is known. Practically, 
this is done by running preliminary analyses on coarse meshes to obtain an 
approximative crack path and heuristically refining the mesh around this path 
to increase accuracy. 

Typically, this refinement does not rely on sound error measures, thus the 
heuristically chosen mesh is in general inadequately suited and can cause large 
inaccuracies in the crack growth path, especially around holes, which can lead 
to non-conservative estimates of the safe life of the structure. 

Thus, it is clear that although XFEM simplifies the treatment of cracks 
compared to the standard FEM by lifting the burden of a geometry conforming 
mesh, it still requires iterations and associated user intervention and it employs 
heuristics which are detrimental to the robustness and accuracy of the simulation 
process. 

It would be desirable to minimise the changes to the mesh topology, thus 
user intervention, while ensuring the stress fields and SIFs are accurately com- 
puted at each crack growth step. This paper is one more step in attempting to 
control the discretization error committed by enriched FE approximations, de- 
crease human intervention in damage tolerance assessment of complex industrial 
structures and enhance confidence in the results by providing enriched FEMs 
with sound error estimators which will guarantee a predetermined accuracy level 
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and suppress recourse to manual iterations and heuristics. 

The error assessment procedures used in finite element analysis are well 



known and can be usually classified into different families ( Ainsworth and Oden 
2000): residual based error estimators, recovery based error estimators, dual 
techniques, etc. Residual based error estimators were substantially improved 
with the introduction of the residual equilibration by [Ainsworth and Oden| 
(2000). These error estimators have a strong mathematical basis and have been 



frequently used to obtain lower and upper bounds of the error (Ainsworth and 



Oden 2000 Diez et al. 2004). Recovery based error estimates were first intro- 



duced by Zienkiewicz and Zhu ( 1987 ) and are often preferred by practitioners 
because they are robust and simple to use and provide an enhanced solution. 
Further improvements were made with the introduction of new recovery pro- 
cesses such as the superconvergent patch recovery (SPR) technique proposed 
by Zienkiewicz and Zhu ( 1992a|b ). Dual techniques based on the evaluation of 
two different fields, one compatible in displacements and another equilibrated 
in stresses, have also been used by Pereira et al. ( 1999 1 to obtain bounds of the 
error. Herein we are going to focus on recovery based techniques which follow 



the ideas of the Zienkiewicz-Zhu (ZZ) error estimator proposed by Zienkiewicz 



and Zhu (1987). 



The literature on error estimation techniques for mesh based partition of 
unity methods is still scarce. One of the first steps in that direction was made in 



the context of the Generalized Finite Element Method (GFEM) by Strouboulis 



et al. (2001 1. The authors proposed a recovery-based error estimator which pro- 
vides good results for /i-adapted meshes. In a later work two new a posteriori 
error estimators for GFEM were presented (Strouboulis et al. 2006). The first 
one was based on patch residual indicators and provided an accurate theoret- 
ical upper bound estimate, but its computed version severely underestimated 
the exact error. The second one was an error estimator based on a recovered 
displacement field and its performance was closely related to the quality of the 
GFEM solution. This recovery technique constructs a recovered solution on 
patches using a basis enriched with handbook functions. 

In order to obtain a recovered stress field that improves the accuracy of 
the stresses obtained by XFEM, Xiao and Karihaloo (2006) proposed a moving 
least squares fitting adapted to the XFEM framework which considers the use 
of statically admissible basis functions. Nevertheless, the recovered stress field 
was not used to obtain an error indicator. 



Pannachet et al. (2009) worked on error estimation for mesh adaptivity in 



XFEM for cohesive crack problems. Two error estimates were used, one based 
on the error in the energy norm and another that considers the error in a local 
quantity of interest. The error estimation was based on solving a series of local 
problems with prescribed homogeneous boundary conditions. 



Panetier et al. (2010) presented an extension to enriched approximations 



of the constitutive relation error (CRE) technique already available to evaluate 
error bounds in FEM ( Ladeveze and Peile| [2005 1 . This procedure has been used 
to obtain local error bounds on quantities of interest for XFEM problems. 



Bordas and Duflot ( 2007 1 and Bordas et al. ( 2008 1 proposed a recovery based 



error estimator for 2D and 3D XFEM approximations for cracks known as the 
extended moving least squares (XMLS). This method intrinsically enriched an 
MLS formulation to include information about the singular fields near the crack 
tip. Additionally, it used a diffraction method to introduce the discontinuity in 
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the recovered field. This error estimator provided accurate results with efFec- 
tivity indices close to unity (optimal value) for 2D and 3D fracture mechanics 
problems. Later, Duflot and Bordas (2008 1 proposed a global derivative recovery 



formulation extended to XFEM problems. The recovered solution was sought 
in a space spanned by the near tip strain fields obtained from differentiating 
the Westergaard asymptotic expansion. Although the results provided by this 



technique were not as accurate as those in Bordas and Duflot (2007), Bordas 



et al. (20081, they were deemed by the authors to require less computational 



power. 



Rodenas et al. ( 2008 ) presented a modification of the superconvergent patch 



recovery (SPR) technique tailored to the XFEM framework called SPRxfem- 
This technique was based on three key ingredients: (1) the use of a singu- 
lar+smooth stress field decomposition procedure around the crack tip similar 



to that described by Rodenas et al. (2006 1 for FEM; (2) direct calculation of re- 
covered stresses at integration points using the partition of unity, and (3) use of 
different stress interpolation polynomials at each side of the crack when a crack 
intersects a patch. In order to obtain an equilibrated smooth recovered stress 
field, a simplified version of the SPR-C technique presented by [Rodenas et al.| 
( 2007 ) was used. This simplified SPR-C imposed the fulfilment of the boundary 



equilibrium equation at boundary nodes but did not impose the satisfaction of 
the internal equilibrium equations. 



The numerical results presented by Bordas and Duflot (2007), Bordas et al. 



( 2008 ) and Rodenas et al. ( 2008 1 showed promising accuracy for both the XMLS 



and the SPRxfem techniques. It is apparent in those papers that SPRxfem led 
to effectivity indices remarkably close to unity. Yet, these papers do not shed 
any significant light on the respective roles played by two important ingredients 
in those error estimators, namely: 

1. The statistical admissibility of the recovered solution; 

2. The enrichment of the recovered solution. 

Moreover, the differences in the test cases analysed and in the quality mea- 
sures considered in each of the papers makes it difficult to objectively compare 
the merits of both methods. 

The aim of this paper is to assess the role of statistical admissibility and en- 
richment of the recovered solution in recovery based error estimation of enriched 
finite element approximation for linear elastic fracture. To do so, we perform 
a systematic study of the results obtained when considering the different fea- 
tures of two error estimation techniques: XMLS and SPR-CX. SPR-CX is an 



enhanced version of the SPRxfem presented by Rodenas et al.| (|2008| for 2D 
problems. Advantages and disadvantages of each method are also provided. 

The outline of the paper is as follows. In Section 2, the XFEM is brieffy 
presented. Section 3 deals with error estimation and quality assessment of the 
solution. In Sections 4 and 5 we introduce the error estimators used in XFEM 
approximations: the SPR-CX and the XMLS, respectively. Some numerical 
examples analysing both techniques and the effect of the enrichment functions 
in the recovery process are presented in Section 6. Finally some concluding 
remarks are provided in Section 7. 
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2 Reference problem and XFEM solution 



Let us consider a 2D linear elastic fracture mechanics (LEFM) problem on a 
bounded domain 51 C M^. The unknown displacement field u is the solution of 
the boundary value problem 



V • tT(u) + b = 
a(u) 11 = 1 
ct(u) n = 
u = u 



in n 

on Fjv 
on Tc 

on r_D 



(1) 
(2) 
(3) 
(4) 



where F^v and T £> are the Neumann and Dirichlet boundaries and Tq represents 
the crack faces, with = Fat U Fd U Tc and F^r n F^ H Tc — 0. h are the 
body forces per unit volume, t the tractions applied on Tn (being n the normal 
vector to the boundary) and u the prescribed displacements on F/j. The weak 
form of the problem reads: Find u £ V such that: 



yveV a(u,v) ?(v). 



(5) 



where V is the standard test space for elasticity problems such that = {v | v G 
v|rj3 (x) = 0, V discontinuous on F^}, and 



a(u, v) := / ct(u) : e(v)df2 = / o-(u) : S : ct(v' 
Jn Jn 

:= / b • vdn + t • vdF, 



(6) 
(7) 



where S is the compliance tensor, cr and e represent the stress and strain oper- 
ators. 

LEFM problems are denoted by the singularity present at the crack tip. The 
following expressions represent the first term of the asymptotic expansion which 
describes the displacements and stresses for combined loading modes I and II 



in 2D. These expressions can be found in the literature (Szabo and Babuska 



1991 Rodenas et al. 2008 1 and are reproduced here for completeness: 
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where r and are the crack tip polar coordinates, Ki and Ku are the stress 
intensity factors for modes I and II, fj, is the shear modulus, and k the Kolosov's 
constant, defined in terms of the parameters of material E (Young's modulus) 
and V (Poisson's ratio), according to the expressions: 



E 



3 — 4u plane strain 



^yL + u plane stress 

This type of problem is difficult to model using a standard FEM approxima- 
tion as the mesh needs to explicitly conform to the crack geometry. With the 
XFEM the discontinuity of the displacement field along the crack faces is intro- 
duced by adding degrees of freedom to the nodes of the elements intersected by 
the crack. This tackles the problem of adjusting the mesh to the geometry of 



the crack (Moes et al. 19991 Stolarska et al. 2001). Additionally, to describe 



the solution around the crack tip, the numerical model introduces a basis that 
spans the near tip asymptotic fields. The following expression is generally used 
to interpolate the displacements at a point of coordinates x accounting for the 
presence of a crack tip in a 2D XFEM approximation: 



u,(x) = + E N,{^)H{^)h, + N,^i^) Y (10) 



where Ni are the shape functions associated with node i, represent the con- 
ventional nodal degrees of freedom, hj are the coefficients associated with the 
discontinuous enrichment functions, and those associated with the functions 
spanning the asymptotic field. In the above equation, I is the set of all the nodes 
in the mesh, Ai is the subset of nodes enriched with crack tip functions, and J' 
is the subset of nodes enriched with the discontinuous enrichment (see Figure [T|). 
In (10 1, the Heaviside function H, with unitary modulus and a change of sign 



on the crack face, describes the displacement discontinuity if the finite element 
is intersected by the crack. The Eg are the set of branch functions used to rep- 
resent the asymptotic expansion of the displacement field around the crack tip 
seen in ([8|. For the 2D case, the following functions are used (Belytschko and 



Black 19991: 



{Ei{r,(j))} = V?'|sin-,cos-,sin-sin0,cos-sin0| (11) 
The main features of the XFEM implementation considered to evaluate the 



numerical results is described in detail in Rodenas et al. ( 2008 1 and can be 
summarized as follows: 

• Use of bilinear quadrilaterals. 

• Decomposition of elements intersected by the crack into integration sub- 
domains that do not contain the crack. Alternatives which do not required 



this subdivision are proposed by Ventura (20061, Natarajan et al. (20101 



• Use of a quasi-polar integration with a 5 x 5 quadrature rule in triangular 
subdomains for elements containing the crack tip. 
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□ F, enrichment 
O enrichment 

A H+F^ em'ichment 

— crack 



Figure 1: Classification of nodes in XFEM. Fixed enrichment area of radius 



No correction for blending elements. Methods to address blending errors 



Tarancon et al. (20091. 



are proposed by Chessa et al. (20031, Grade et al. (2008), Fries (20081, 



2.1 Evaluation of stress intensity factors 

The stress intensity factors (SIFs) in LEFM represent the amplitude of the sin- 
gular stress fields and are key quantities of interest to simulate crack growth in 
LEFM. Several post-processing methods, following local or global (energy) ap- 
proaches, are commonly used to extract SIFs (Banks-Sills 19911 or to calculate 
the energy release rate G. Energy or global methods are considered to be the 
most accurate and efficient methods ( Banks-Sillsl 19911 Li et al. 1985). Global 



methods based on the equivalent domain integral (EDI methods) are specially 
well-suited for FEM and XFEM analyses as they are easy to implement and can 
use information far from the singularity. In this paper, the interaction integral 



as described in Shih and Asaro ( 1988 1, Yau et al. ( 1980 1 has been used to extract 
the SIFs. This technique provides Ki and Ku for problems under mixed mode 
loading conditions using auxiliary fields. Details on the implementation of the 



interaction integral can be found for example in Moes et al. (1999), Rodenas 



et al. (2008), Giner et al. (2005). 



3 Error estimation in the energy norm 

The approximate nature of the FEM and XFEM approximations implies a dis- 
cretization error which can be quantified using the error in the energy norm for 
the solution ||e|| = ||u — u'^jj. To obtain an estimate of the discretization error 
||ees||, in the context of elasticity problems solved using the FEM, the expression 
for the ZZ estimator is defined as (Zienkiewicz and Zhu 1987[ ) (in matrix form): 



/ {a* - a'^f B-' {cT* - a"") dCi 
Jn 



(12) 
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or alternatively for the strains 

Jn 

where the domain 51 refers to the full domain of the problem or a local sub- 
domain (element), cr'* represents the stress field provided by the FEM, cr* is 
the recovered stress field, which is a better approximation to the exact solution 
than (T^ and D is the elasticity matrix of the constitutive relation cr = De. 

The recovered stress field cr* is usually interpolated in each element using 
the shape functions N of the underlying FE approximation and the values of 
the recovered stress field calculated at the nodes a* 



a*(x) = ^7V,(x)^:(x,), (13) 

1=1 

where is the number of nodes in the element under consideration and o'*(xi) 
are the stresses provided by the least squares technique at node i. The com- 
ponents of a* are obtained using a polynomial expansion, aij = pa (with 
j — xx,yy,xy), defined over a set of contiguous elements connected to node i 
called patch, where p is the polynomial basis and a are the unknown coefficients. 

The ZZ error estimator is asymptotically exact if the recovered solution used 
in the error estimation converges at a higher rate than the finite element solution 



(Zienkiewicz and Zhu 1992b). 

In this paper, we are interested in the role played by statical admissibility 
and enrichment of the recovered solution in estimating the error committed 
by XFEM. To do so, we study the performance of two recovery-based error 
estimators, which exhibit different features, that have been recently developed 
for XFEM: 

• The SPR-CX derived from the error estimator developed by |R6denas at al.] 
(j2008j) and summarized in Section |4] (two other versions of the SPR-CX 
technique have been also considered); 



The XMLS proposed by Bordas and Duflot (20071 and summarized in 
Section [s] 



Both estimators provide a recovered stress field in order to evaluate the 



estimated error in the energy norm by means of the expression shown in ( 12 I 



4 SPR-CX error estimator 



The SPR-CX error estimator is an enhancement of the error estimator first 
introduced by Rodenas et al. (20081, which incorporates the ideas proposed 



Rodenas et al. ( 2007 1 to guarantee the exact satisfaction of the equilibrium 



locally on patches. In Rodenas et al. (20081 a set of key ideas are proposed to 



modify the standard SPR by Zienkiewicz and Zhu (1992a I, allowing its use for 
singular problems. 

The recovered stresses cr* are directly evaluated at an integration point x 
through the use of a partition of unity procedure, properly weighting the stress 
interpolation polynomials obtained from the different patches formed at the 
vertex nodes of the element containing x: 
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a*(x)=^7V,(xK(x), (14) 

1=1 

where Ni are the shape functions associated to the vertex nodes n„. 

One major modification is the introduction of a spHtting procedure to per- 
form the recovery. For singular problems the exact stress field cr is decomposed 
into two stress fields, a smooth field ct^^q and a singular field (Tsing- 

— ^ smo ~t- (^sing- (-^^) 

Then, the recovered field cr* required to compute the error estimate given 



in (12) can be expressed as the contribution of two recovered stress fields, one 
smooth cr*„jQ and one singular cr*ging'- 

= (^*smo + ^ling- (16) 

For the recovery of the singular part, the expressions in (|9| which describe 
the asymptotic fields near the crack tip are used. To evaluate cr*^^^ from (j9]) 
we first obtain estimated values of the stress intensity factors Ki and Kn using 



the interaction integral as indicated in Section |2.1[ The recovered part cr^^j^g is 



an equilibrated field as it satisfies the internal equilibrium equations 
Once the field cr^i^g has been evaluated, an FE approximation to t 
part cr^yno can be obtained subtracting cr*^,^^ from the raw FE field: 



Then, the field cr*^^, is evaluated applying an SPR-C recovery procedure 
over the field cr'^rno- 

For patches intersected by the crack, the recovery technique uses different 
stress interpolation polynomials on each side of the crack. This way it can 
represent the discontinuity of the recovered stress field along the crack faces, 
which is not the case for SPR that smoothes out the discontinuity (see, e.g. 
Bordas and Duflot| ( p007| ) 



In order to obtain an equilibrated recovered stress field cr*^^^, the SPR-CX 
enforces the fulfilment of the equilibrium equations locally on each patch. The 
constraint equations are introduced via Lagrange multipliers into the linear 
system used to solve for the coefficients of the polynomial expansion of the 
recovered stresses on each patch. These include the satisfaction of the: 

• Internal equilibrium equations. 

• Boundary equilibrium equation: A point collocation approach is used to 
impose the satisfaction of a second order approximation to the tractions 
along the Neumann boundary. 

• Compatibility equation: This additional constraint is also imposed to fur- 
ther increase the accuracy of the recovered stress field. 

To evaluate the recovered field, quadratic polynomials have been used in 
the patches along the boundary and crack faces, and linear polynomials for the 
remaining patches. As more information about the solution is available along 
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the boundary, polynomials one degree higher are useful to improve the quality 
of the recovered stress field. 

The enforcement of equilibrium equations provides an equilibrated recovered 
stress field locally on patches. However, the process used to obtain a continuous 
field a* shown in ( 14 1 introduces a small lack of equilibrium as explained in 



Rodenas et al. (2010al. The reader is referred to Rodenas at al. (2010al, Diez 



et al. (2007) for details. 



5 XMLS error estimator 



In the XMLS the solution is recovered through the use of the moving least 
squares (MLS) technique, developed by mathematicians to build and fit sur- 
faces. The XMLS technique extends the work oflTabbara et al.|(|1994|) for FEM 



to enriched approximations. The general idea of the XMLS is to use the dis- 



placement solution provided by XFEM to obtain a recovered strain field ( Bordas 



and Duflot 2007 Bordas et al. 20081. The smoothed strains are recovered from 



the derivative of the MLS-smoothed XFEM displacement field: 

u*(x)= ^ V^,(x)u,f 

£*(x)^ ^ V.(^,(x)u,''), 



(18) 
(19) 



where the are the raw nodal XFEM displacements and V'j(x) are the MLS 
shape function values associated with a node i at a point x. A/'x is the set of 
nodes in the domain of infiuence of point x, Vs is the symmetric gradient 
operator, u* and e* are the recovered displacement and strain fields respectively. 
At each point x^ the MLS shape functions "0^ are evaluated using weighting 
functions uji and an enriched basis p(xi). The total Ux non-zero MLS shape 
functions at point x are evaluated as: 

(x)p(x))^p(x,)tj,(x) 



where A(x) 



(V',(x))i<,<„^ ^ (A-^(x)p(x))^ p(x,)t^.(x) (20) 
S"=i '^i(x)p(^i)p(^j)"^ is 3, matrix to be inverted at every point 



X (see [Bordas and Dufiot] ( |2007[ ), [Bordas et all ( |2008[ ) for further details). 

For each supporting point x^, the weighting function uj^ is defined such that: 



UJl{s) = fi{s) 



6s2 



if \s\ < 1 
if Isl > 1 



(21) 



where s is the normalized distance between the supporting point x^ and a point 
X in the computational domain. In order to describe the discontinuity, the 
distance s in the weight function defined for each supporting point is modified 
using the diffraction criterion (Belytschko et al. 19961. The basic idea of this 



criterion is depicted in Figure [2] The weight function is continuous except across 
the crack faces since the points at the other side of the crack are not considered 
as part of the support. Near the crack tip the weight of a node i over a point 
of coordinates x diminishes as the crack hides the point. When the point x is 
hidden by the crack the following expression is used: 



di 



(22) 



10 



where di is the radius of the support. 




Figure 2: Diffraction criteria to introduce the discontinuity in the XMLS ap- 
proximation. 

The MLS shape functions can reproduce any function in their basis. The 
basis p used is a Hnear basis enriched with the functions that describe the first 
order asymptotic expansion at the crack tip as indicated in ( |11[ ) : 

p=[l,x,y, [F,{r, <^), F2ir, c^),F^{r, c^),Fi{r, 0)]] (23) 

Note that although the enriched basis can reproduce the singular behaviour 
of the solution around the crack tip, the resulting recovered field not necessarily 
satisfies the equilibrium equations. 



6 Numerical results 

In this section, numerical experiments are performed to verify the behaviour of 
both XFEM recovery based error estimators considered in this paper. 

Babuska et al. ( 1994a|b 19971 proposed a robustness patch test for quality 
assessment of error estimators. However, the use of this test is not within the 
scope of this paper and furthermore, to the authors' knowledge, it has not been 
used in the context of XFEM. Therefore, the more traditional approach of using 
benchmark problems is considered here to analyse the response of the different 
error estimators. 

The accuracy of the error estimators is evaluated both locally and globally. 
This evaluation has been based on the effectivity of the error in the energy norm, 
which is quantified using the effectivity index 9 defined as: 



(24) 



When the enhanced or recovered solution is close to the analytical solution 
the effectivity approaches the theoretical value of 1, which indicates that it is a 
good error estimator, i.e. the approximate error is close to the exact error. 

To assess the quality of the estimator at a local level, the local effectivity 
D, inspired on the robustness index found in Babuska et al. (1994a), is used. 



For each element e, D represents the variation of the effectivity index in this 
element, 9'^ , with respect to the theoretical value (the error estimator can be 
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considered to be of good quality if it yields D values close to zero) . D is defined 
according to the following expression, where superscripts indicate the element 
e: 



D = 61^ - 1 



if 
if 



> 1 

' < 1 



with 



(25) 



To evaluate the overall quality of the error estimator we use the global ef- 
fectivity index 9, the mean value mdZ^I) and the standard deviation <J {D) of 
the local effectivity. A good quality error estimator yields values of 9 close to 
one and values of m {\D\) and a (D) close to zero. 

The techniques can be used in practical applications. However, in order to 
properly compare their performance we have used an academic problem with 
exact solution. In the analysis we solve the Westergaard problem (Gdoutos 



1993) as it is one of the few problems in LEFM under mixed mode with an 



analytical solution. In |Giner et al. (20051, Rodenas et al. (20081 we can find 
explicit expressions for the stress fields in terms of the spatial coordinates. In the 
next subsection we show a description of the Westergaard problem and XFEM 



model, taken from Rodenas et al. ( 2008 1 and reproduced here for completeness. 



6.1 Westergaard problem and XFEM model 

The Westergaard problem corresponds to an infinite plate loaded at infinity 
with biaxial tractions axoo = o'yoo = <^oo and shear traction Too, presenting a 
crack of length 2a as shown in Figure [3j Combining the externally applied loads 
we can obtain different loading conditions: pure mode I, II or mixed mode. 



a. 
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-I 
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Figure 3: Westergaard problem. Infinite plate with a crack of length 2a un- 
der uniform tractions CToo (biaxial) and Too- Finite portion of the domain CIq, 
modelled with FE. 

In the numerical model only a finite portion of the domain (a — 1 and 
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6 = 4 in Figure [3| is considered. The projection of the stress distribution 
corresponding to the analytical Westergaard solution for modes I and II, given 
by the expressions below, is applied to its boundary: 
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(27) 

where the stress fields are expressed as a function of x and y, with origin at the 
centre of the crack. The parameters t, m, n and cj) are defined as 



t = {x + lyY — = [x^ ~ ~ o?') + ii'^xy) = m + in 
m = Re{t) = Re(z2 - a^) = - - 
n = Im(t) = (z^ - a^) = 2a;j/ 

= Arg(f) = Arg(m — in) with (/) G tt, tt] , = —1 
For the problem analysed, the exact value of the SIF is given as 



(28) 
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(29) 



Three different loading configurations corresponding to the pure mode I 
(coo = 100, Too — 0), pure mode II (uoo = 0, Too = 100), and mixed mode 
(coo = 30, Too — 90) cases of the Westergaard problem are considered. The 
geometric models and boundary conditions are shown in Figure |4] 

Bilinear elements are used in the models, with a singular + smooth decompo- 
sition area of radius p = 0.5 equal to the radius re used for the fixed enrichment 
area (note that the splitting radius p should be greater or equal to the enrich- 
ment radius re, Rodenas et al. ( 2008[ )). The radius for the Plateau function 



used to extract the SIF is = 0.9. Young's modulus is i? = 10^ and Poisson's 
ratio V — 0.333. For the numerical integration of standard elements we use a 
2x2 Gaussian quadrature rule. For split elements we use 7 Gauss points in 



each triangular subdomain, and a 5 x 5 quasipolar integration (Bechet et al. 



2005 ) in the subdomains of the element containing the crack tip. 
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Regarding global error estimation, the evolution of global parameters in 
sequences of uniformly refined structured (Figurejs]) and unstructured (Figure[6]) 
meshes is studied. In the first case, the mesh sequence is defined so that the crack 
tip always coincided with a node. For a more general scenario, this condition 
has not been applied for the unstructured meshes. 




Mesh 1 Mesh 2 Mesh 3 Mesh 4 Mesh 5 

288 elem. 512elem. 800 elem. 3200 clem. 12800 elem. 



Figure 5: Sequence of structured meshes. 



6.2 Mode I and structured meshes 

The first set of results presented in this study are for the Westergaard problem 
under mode I load conditions. Figure [7] shows the values for the local effectivity 
index D for the first mesh in the sequence analysed. In the figure, the size of 
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Meshl 
273 elem. 



Mesh 2 
525 clem. 



Mesh 3 
818 elem. 



Mesh 4 
3272 elem. 



Mesh 5 
13088 dcm. 



Figure 6: Sequence of unstructured meshes. 



the enrichment area is denoted by a circle. It can be observed that far from 
the enrichment area both techniques yield similar results, however, the XMLS 
technique exhibits a higher overestimation of the error close to the singularity. 



I 



I 



a. SPR-CX 



b. XMLS 



Figure 7: Mode I, structured mesh 1. Local effectivity index D (ideal value 
D = 0). The circles denote the enriched zones around the crack tips. The 
superiority of the statically admissible recovery (SPR-CX) compared to the 
standard XMLS is clear. 

The same behaviour is also observed for more refined meshes in the sequence. 
Figure [8] shows a zoom in the enriched area of a finer mesh where, as before, a 
higher overestimation of the error can be observed for the results obtained with 
the XMLS error estimator. In this example and in order to study the evolution 
of the accuracy of the recovered stress field when considering different features 
in the recovery process, we consider two additional versions: 

• The first considers a recovery procedure which enforces both internal and 
boundary equilibrium, but does not include the singular+smooth splitting 
technique (SPR-C). This approach is very similar in form to conventional 
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SPR-based recovery techniques widely used in FEM. It is well-known that 
this type of recovery process will produce unreliable results when the stress 
field contains a singularity because the polynomial representation of the 
recovered stresses is not able to describe the singular field (e.g. Bordas 



and Duflot (2007)). 



The second version of the error estimator performs the splitting but does 
not equilibrate the recovered field (SPR-X). As previously commented, 
""sing equilibrated field but (t*,^^ is not equilibrated. The aim of 

this second version is to assess the infiuence of enforcing the equilibrium 
constraints. 



Table 1 summarizes the main features of the different recovery procedures we 
considered. 





Singular 


Equilibrated sin- 


Locally equili- 




functions 


gular functions 


brated field 


SPR-X 


yes 


yes 


no 


SPR-C 


no 


not applicable 


yes 


SPR-CX 


yes 


yes 


yes 


XMLS 


yes 


no 


no 



Table 1 : Comparison of features for the different recovery procedures 



The performance of SPR-X in Figure [8]is poor, especially along the Neumann 
boundaries, where the equilibrium equations are not enforced, as they are in the 
SPR-CX. Particularly interesting are the results for the SPR-C technique, where 
considerable overestimation of the exact error is observed in the whole enriched 
region. This is due to the inability of polynomial functions to reproduce the 
near-tip fields, even relatively far from the crack tip. 

The above compared the relative performance of the methods locally. It is 
also useful to have a measure of the global accuracy of the different error estima- 
tors. Figure [gjleft) shows the convergence of the estimated error in the energy 
norm ||ees|| (see equation ( [l2| ) evaluated with the proposed techniques, the 
convergence of the exact error ||e|[^is shown for comparison. It can be observed 
that SPR-CX leads to the best approximation to the exact error in the energy 
norm. The best approach is that which captures as closely as possible the exact 
error. Non-singular formulations such as SPR-C greatly overestimate the error, 
especially close to the singular point. Only SPR-C suffers from the characteris- 
tic suboptimal convergence rate (factor 1/2, associated with the strength of the 
singularity) 

In order to evaluate the quality of the recovered field cr* obtained with each 
of the recovery techniques, the convergence in energy norm of the approximate 
error on the recovered field ||e*|| = ||u — u*|| is compared for the error estimators 
in Figure [gjright). It can be seen that the XMLS, the SPR-CX and the SPR-X 

^The exact error measures the error of the raw XFEM compared to the exact solution. The 
theoretical convergence rate is 0{hP), h being the element size and p the polynomial degree 
(1 for linear elements,...), or (!?(— dof if we consider the number of degrees of freedom. 
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a. SPR-CX 



b. XMLS 





c. SPR-X 



d. SPR-C 



Figure 8: Mode I. structured mesh 4. Local efTectivity index D (ideal value 
D = 0). The circles denote the enriched zones around the crack tips. The 
results show the need for the inclusion of the near-tip fields in the recovery 
process (SPR-X is far superior to SPR-C). It is also clear that enforcing statical 
admissibility (SPR-CX) greatly improves the accuracy along the crack faces and 
the Neumann boundaries. Also notice that XMLS leads to slightly improved 
effectivities compared to SPR-CX around the crack faces, between the enriched 
zone and the left boundary. SPR-C is clearly not able to reproduce the singular 
fields, which is shown by large values of D inside the enriched region. 



provide convergence rates higher than the convergence rate for the exact error 
||e||, which is an indicator of the asymptotic exactness of the error estimators 
based on these recovery techniques ( [Zienkiewicz and Zhu, ,1992b| ) . 

Figure 10 shows the evolution of the parameters 9, m{\D\) y a (D) with 
respect to the number of degrees of freedom for the error estimators. It can 
be verified that although the SPR-C technique includes the satisfaction of the 
equilibrium equations, the effectivity of the error estimator does not converge 
to the theoretical value {6 = 1). This is due to the absence of the near-tip 
fields in the recovered solution. The influence of introducing singular functions 
in the recovery process can be observed in the results provided by SPR-X. In 
this case, the convergence is obtained both locally and globally. This shows 
the importance of introducing a description of the singular field in the recovery 
process, i.e. of using extended recovery techniques in XFEM. 

Regarding the enforcement of the equilibrium equations, we can see that 
an equilibrated formulation (SPR-CX) leads to better effectivities compared to 
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Figure 9: Mode I, structured meshes: Convergence of the estimated error in 
the energy norm Heg^H (left). All methods which include the near-tip fields 
do converge with the optimal convergence rate of 1/2. The best method is 
SPR-CX. XMLS performs worse than SPR-X for coarse meshes, but it becomes 
equivalent to SPR-CX for fine meshes whilst the SPR-X error increases slightly. 
Convergence of the error for the recovered field ||e*g|| (right). For the recovered 
field, SPR-CX is the best method closely followed by SPR-X. XMLS errors 
are half an order of magnitude larger than SPR-CX/SPR-X, but superior to the 
non-enriched recovery method SPR-C. For all methods, except SPR-C, the error 
on the recovered solution converges faster than the error on the raw solution, 
which shows that the estimators are asymptotically exact. 



other non-equilibrated configurations (SPR-X, XMLS), which is an indication 
of the advantages associated with equilibrated recoveries in this context. 

Still in Figure [TOj the results for the SPR-CX and XMLS show that the 
values of the global effectivity index are close to (and less than) unity for the 
SPR-CX estimator, whilst the XMLS is a lot less effective and shows effectivities 
larger than 1. 

The next step in our analysis is the verification of the convergence of the 
mean value and standard deviation of the local effectivity index D. From a 
practical perspective, one would want the local effectivity to be equally good, 
on average, everywhere inside the domain; one would also want this property 
to improve with mesh refinement. In other words, the average local effectivity 
m and its standard deviation a should decrease with mesh refinement: m and 
cr measure the average behaviour of the method and how far the results deviate 
from the mean, i.e. how spatially consistent they are. The idea is to spot areas 
where there may be compensation between overestimated and underestimated 
areas that could produce an apparently 'accurate' global estimation. It can be 
confirmed that, although the curves for m{\D\) and (t{D) for both estimators 
tend to zero with mesh refinement, the SPR-CX is clearly superior to the XMLS. 
These results can be explained by the fact that the SPR-CX enforces the ful- 
filment of the equilibrium equations in patches, and evaluates the singular part 
of the recovered field using the equilibrated expression that represents the first 
term of the asymptotic expansion, whereas the XMLS enrichment uses only a 
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Figure 10: Global indicators 0, m(|D|) and (J (D) for mode I and structured 
meshes. 



set of singular functions that would be able to reproduce this term but are not 
necessarily equilibrated. 

Although the results for the SPR-CX proved to be more accurate, the XMLS- 
type techniques will prove useful to obtain error bounds. There is an increasing 
interest in evaluating upper and lower bounds of the error for XFEM approxima- 
tions. Some work has been already done to obtain upper bounds using recovery 
techniques as indicated in Rodenas et al. (2010a), where the authors showed 



that an upper bound can be obtained if the recovered field is continuous and 
equilibrated. However, they demonstrated that due to the use of the conjoint 
polynomials process to enforce the continuity of cr*, some residuals in the equi- 
librium equations appear, and the recovered field is only nearly equilibrated. 
Then, correction terms have to be evaluated to obtain the upper bound. An 
equilibrated version of the XMLS could provide a recovered stress field which 
would be continuous and equilibrated, thus facilitating the evaluation of the up- 
per bound. Initial results for this class of techniques can be found in |R6denas] 



et al. (2009 20 10b I. 



In this first example we have shown the effect of using the SPR-C and SPR-X 
recovery techniques in the error estimation. Considering that these two tech- 
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niques are special cases of the SPR-CX with inferior results, for further examples 
we will focus only on the SPR-CX and XMLS techniques. 



6.3 Mode II and structured meshes 

Figure [TT] presents the results considering mode II loading conditions for the 
local effectivity index D on the fourth mesh of the sequence. Similarly to the 
results for mode I, the XMLS estimator presents a higher overestimation of the 
error near the enrichment area. This same behaviour is observed for the whole 
set of structured meshes analysed under pure mode II. 




a. SPR-CX b. XMLS 



Figure 11: Mode II, structured mesh 4. Local effectivity index D (ideal value 
D = Q). SPR-CX leads to much better local effectivities than XMLS. It is also 
remarkable that the effectivities are clearly worse in the enriched region (circle) 
for both the XMLS and the SPR-CX, and that this effect is more pronounced 
in the former method. Note the slightly worse results obtained by SPR-CX 
around the crack faces between the boundary of the enriched region and the left 
boundary, as in mode I. 



As for the mode I case, the evolution of global accuracy parameters in mode 
II exhibits the same behaviour seen for mode I loading conditions. Figure 
shows the results for the convergence of the error in the energy norm. Figure 



shows the evolution of global indicators 6, m {\D\) and a (D). Again, SPR-CX 
provides the best results for the Westergaard problem under pure mode II, using 
structured meshes. 



6.4 Mixed mode and unstructured meshes 

Considering a more general problem. Figure [14] shows the local effectivity index 
D for the fourth mesh in a sequence of unstructured meshes, having a number 
of degrees of freedom (dof ) similar to the mesh represented previously for load 
modes I and II. The XMLS recovered field presents higher overestimation of the 
error around the crack tip which corroborates the results found in previous load 
cases. 
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Figure 12: Mode II, structured meshes: Convergence of the estimated error in 
the energy norm Heg^H (left). The convergence rates are very similar to those 
obtained in the mode I case, and, SPR-CX is still the most accurate method, i.e. 
that for which the estimated error is closest to the exact error. Convergence of 
the error in the recovered field ||e*^|| (right). It can be noticed that the SPR-CX 
error still converges faster than the exact error, but with a lower convergence 
rate (0.59 versus 0.73) than in the mode I case. There is no such difference 
for the XMLS results, which converge at practically the same rate (0.71 versus 
0.72). The error level difference between SPR-CX and XMLS is about half an 
order of magnitude, as in the mode I case. 



Figures [15] and [16] represent the evolution of global parameters for unstruc- 
tured meshes and mixed mode. Once more, the best results for the error esti- 
mates are obtained using the SPR-CX technique. 



7 Conclusions 



The aim of this paper was to assess the accuracy gains provided by 

1. Including relevant enrichment functions in the recovery process; 

2. Enforcing statical admissibility of the recovered solution. 

We focused on two recovery-based error estimators already available for 
LEFM problems using the XFEM. The first technique called SPR-CX is an 
enhancement of the SPR-based error estimator presented by [Rodenas et al.| 
( 2008 1 , where the stress field is split into two parts (singular and smooth) and 



equilibrium equations are enforced locally on patches. The second technique is 
the XMLS proposed by jBordas and Duflot] ( |2007[ ), [Bordas et e^ ( |2008[ ) which 
enriches the basis of MLS shape functions, and uses a diffraction criterion, in 
order to capture the discontinuity along the crack faces and the singularity at 
the crack tip. 

To analyse the behaviour of the ZZ error estimator using both techniques 
and to assess the quality of the recovered stress field, we have evaluated the 
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Figure 13: Global indicators 9, m{\D\) and cr (D) for mode II and structured 
meshes. The results are qualitatively and quantitatively similar to those ob- 
tained in mode I. 
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a. SPR-CX 



b. XMLS 



Figure 14: Mixed mode, unstructured mesh 4. Local effectivity index D (ideal 
value D = 0). Note the improved results compared to the structured case with 
D values ranging from -1.5 to 1.5 as opposed to -4 to 4. 



effectivity index considering problems with an exact solution. Convergence of 
the estimated error in the energy norm and other local error indicators are also 
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Figure 15: Mixed mode, unstructured meshes: Convergence of the estimated 
error in the energy norm ||ees|| (left). Convergence of the error for the recovered 
field ||e*j,|| (right). The results are almost identical to the structured mesh case, 
except for the faster convergence obtained for SPR-CX in the unstructured 
compared to the structured case (0.77 versus 0.59), keeping in mind that optimal 
convergence rates are only formally obtained for structured meshes. 



evaluated. To analyse the infiuence of the special features introduced in the 
recovery process we have also considered two additional configurations: SPR-C 
and SPR-X. 

The results indicate that both techniques, SPR-CX and XMLS, provide error 
estimates that converge to the exact value and can be considered as asymptot- 
ically exact. Both techniques could be effectively used to estimate the error in 
XFEM approximations, while other conventional recovery procedures which do 
not include the enrichment functions in the recovery process have proved not to 
converge to the exact error. This shows (albeit only numerically) the need for 
the use of extended recovery techniques for accuracy assessment in the XFEM 
context. 

Better results are systematically obtained when using the SPR-CX to re- 
cover the stress field, specially in the areas close to the singular point. For 
all the different test cases analysed, the XMLS produced higher values of the 
effectivity index in the enriched area, where the SPR-CX proved to be more 
accurate. This can be ascribed to the fact that the SPR-CX technique recovers 
the singular part of the solution using the known equilibrated exact expressions 
for the asymptotic fields around the crack tip and enforces the fulfilment of 
the equilibrium equations on patches. Further work currently in progress will 
include the development of an equilibrated XMLS formulation, which could pro- 
vide a continuous and globally equilibrated recovered stress field which can then 
be used to obtain upper bounds of the error in energy norm. 

The aim of our project is to tackle practical engineering problems such as 
those presented in Bordas and Moran (20061, Wyart et al. (20071 where the 
accuracy of the stress intensity factor is the target. The superiority of SPR-CX 
may then be particularly relevant. To minimise the error on the stress intensity 
factors we will target this error directly, through goal-oriented error estimates. 
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Figure 16: Global indicators 6, m {\D\) and a (D) for mixed mode and unstruc- 
tured meshes. Notice that the XMLS performs more closely to SPR-CX using 
those indicators, for unstructured than for structured meshes, especially when 
measuring the mean value of the effectivities m{\D\). 

This will be reported in a forthcoming publication. However, although the SPR- 
CX results are superior in general, an enhanced version of the XMLS technique 
presented in this paper where we enforce equilibrium conditions could result 
useful to evaluate upper error bounds when considering goal-oriented error es- 
timates, as it directly produces a continuous equilibrated recovered stress field. 

Our next step will be the comparison of available error estimators in three 
dimensional settings in terms of accuracy versus computational cost to minimise 
the error on the crack path and damage tolerance of the structure. 
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